Two-component few-fermion mixtures in a one-dimensional trap: numerical versus 

analytical approach 
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We explore a few-fermion mixture consisting of two components which are repulsively interacting 
and confined in a one-dimensional harmonic trap. Different scenarios of population imbalance 
ranging from the completely imbalanced case where the physics of a single impurity in the Fermi- 
sea is discussed to the partially imbalanced and equal population configurations are investigated. 
For the numerical calculations the multi-configurational time-dependent Hartree (MCTDH) method 
is employed, extending its application to few-fermion systems. Apart from numerical calculations 
we generalize our Ansatz for a correlated pair wave-function proposed in [J] for bosons to mixtures 
of fermions. From weak to strong coupling between the components the energies, the densities 
and the correlation properties of one-dimensional systems change vastly with an upper limit set 
by fermionization where for infinite repulsion all fermions can be mapped to identical ones. The 
numerical and analytical treatments are in good agreement with respect to the description of this 
crossover. We show that for equal populations each pair of different component atoms splits into 
two single peaks in the density while for partial imbalance additional peaks and plateaus arise for 
very strong interaction strengths. The case of a single impurity atom shows rich behaviour of the 
energy and density as we approach fermionization, and is directly connected to recent experiments 
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I. INTRODUCTION 



Ultracold Fermi gases have been a flourishing research 
field in the last decade after the early experiments on 
40 K [f| and 6 Li [(|. Impressive features like superflu- 
idity according to the Bardeen-Cooper-Schrieffer (BCS) 
theory pj have been realized experimentally Q. Two- 
component ensembles prepared usually in two hypcrfinc 
states or with different atomic species are in the core of 
research for novel phases. Ensembles with imbalanced 
populations js| can then lead to the coexistence of nor- 
mal and superfluid phases [l(| or to the more complex 
phase with a spatially varying gap as proposed by Fulde, 
Ferrell, Larkin and Ovchinnikov (FFLO) [ll|- Highly 
imbalanced Fermi mixtures constitute a very well con- 
trollable system to study polaron physics and recently 
two experiments have achieved measurements of pola- 
ronic behaviour for attractive (l3j and repulsive impurity 
atoms [l4| in goo d ag reement with corresponding theo- 
retical treatments [TJ, |l5[ . An open debate in the field is 
the observation of itinerant ferromagnetism induced by 
interactions [la ]. 

The regime of strong inter-atomic coupling is particu- 
larly interesting, since various spectacular properties of 
highly correlated phases arise. This holds in particu- 
lar for one-dimensional (ID) quantum gases which are 
strongly confined with respect to their transversal modes. 
These quasi-lD tubes are usually equipped with an exter- 
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nal harmonic trap for the longitudinal direction. Apart 
from the well-established toolbox of cold atoms [l7| con- 
taining optical traps and magnetic Feshbach resonances 
to control the trapping geometry and the interaction 
strength respectively, there exist for quasi-lD s etup s also 
the so-called confinement induced resonances [IS] mod- 
ifying the free scattering properties due to the presence 
of the confinement. One of the most peculiar features 
in one-dimension is the mapping between bosons and 
fermions for infinite repulsion which induces an effective 
Pauli exclusion principle to the bosonic species [19]. This 
fermionization of bosons, the so-called Tonks-Girardcau 
(TG) gas has been observed experimentally [20( and more 
recently its highly excited counterpart for attractive in- 
teractions has been probed (2lj . The fermionic physics 
in ID may prove to be even richer in phenomenology as 
the first experiments in confinement-induced molecules 
[22| and spin-imbalanced Fermi gases [23| indicate. On 
theoretical side there have been already several studies 
on the corresponding polaronic physics (24U28| . TG and 
super-TG gases specifically their pairing properties and 
dynamics [29H31| . 

In this work, we investigate a few-fermion system 
composed of two species with the interspecies coupling 
strength being the control parameter which is tuned 
from weak to strong interactions. Methodologically 
we employ two approaches: (a) the numerically exact 
multi-configurational time-dependent Hartree (MCTDH) 
method [321 ] (see Appendix), which has been proven to 
be efficient for few-boson systems (33|, applied here for 
first time on a fermionic few-body ensemble and (b) the 
generalization for mixtures of the explicit wave-function 
based on products of pair-functions which the authors 



2 



introduced for identical bosons in Ref. The latter 

Ansatz has been proven to be efficient for the description 
of the few-body properties in the course of the crossover 
from weak interactions to fermionization, and the exten- 
sion here to few-fermion two-component mixtured is also 
in a good agreement capturing the main behaviour of the 
observables. Our present study is in close connection to 
the experimental work [2| , which has achieved the deter- 
ministic loading of a specific number of fermions into a 
trap with a very high fidelity 0, 13] and accurately con- 
trols the imbalance of the population. The state-of-the 
art of these experiments, allows for a study of specific 
configurations of populations of different hyperfine states 
like those which we aim at here: highly imbalanced with 
one impurity atom surrounded by several (1-4) major- 
ity atoms, a balanced system with equal populations of 
totally each species, and a partially imbalanced ensem- 
ble of eg. 3 majority to 2 minority atoms. For these 
cases we will study the energies and the densities as the 
repulsive interaction strength increases, showing thereby 
the mechanisms and features that characterize these few- 
body ensembles along the crossover to fermionization. 
The case of completely imbalanced population, is already 
measured from the experiment in a work to appear soon 

This article is organized as follows: In Section II we 
present the setup and the explicitly correlated wave func- 
tion which is proposed for its description. In the following 
Section III we analyze and compare the analytical and 
numerical results of the energies for all configurations 
of few-fermion mixtures included in this work. Subse- 
quently we report explicitly on properties of the one-body 
density for the impurity case (Sec. IV) and for the par- 
tially as well as completely balanced systems (Sec. V). 
In the last Section VI we present our conclusions and an 
outlook. The Appendix contains a short description of 
the MCTDH method and remarks on its implementation 
for few-fermion systems. 

II. HAMILTONIAN AND ANSATZ 

The system consists of a two component Fermi gas with 
the two species prepared e.g. in the two lowest hyperfine 
states ( for 6 Li \F = 1/2, m F = ±1/2 > see ref. 0,13]) 
confined in a pure ID harmonic oscillator potential. The 
standard experimental approach is to create quasi-lD 
tubes by using a very strong laser field for the transver- 
sal degrees of freedom compared to the longitudinal one 
[20l [21| . The trap becomes then highly anisotropic with 

the characteristic transversal length scale a± = yj M ^ 

much smaller than the longitudinal one ay = y^^j [ w ± 

(ton) is the transversal (longitudinal) harmonic confine- 
ment frequency]. The typical ratio of the confinement 
frequencies lou/u>j_ < 1/10 is adequate to characterize 
the system as ID. Therefore the transverse degrees of 
freedom are energetically frozen occupying only the cor- 



responding ground states. 

Due to antisymmetry s-wave interactions between 
identical fermions (same hyperfine state) are not possi- 
ble, while p-wave interactions are negligible. Therefore 
we remain with a single interaction parameter, namely 
the interaction strength that between atoms belonging 
to different components. Still this has to take into ac- 
count the quasi-lD nature of the confining potential and 
the resulting confinement induced resonances [181 ] . The 
effective ID interaction strength reads: 

2h 2 a 3D ( icq^i^ V 1 

giD ~ ~W^~ [ ~ V2a ± ) ' 

where cl^d is the 3D s-wave scattering length. When 0313 
which can be tuned via magnetic Feshbach resonances ap- 
proaches the transversal confinement length a±, a con- 
finement induced resonance is encountered. Therefore 
the inter-species interaction strength can be regarded as 
tunable over a very large range. 

With this information at hand we may write the 
rescaled ID harmonic oscillator Hamiltonian: 

1 N d 2 1 N 
jff =-2^^ + 2^ X?+ - 9 £ S( Xi - Xj ) (1) 

i=l 1 i=l i<N M <j<N 

such that lengths and energies are scaled by ay and huj\\ 
respectively, and g = gizj/Swyay. The total number of 
atoms is N and Nm is number of atoms belonging to 
the hyperfine state which has the larger population in 
the trap. We remain therefore with a single parameter g, 
the rescaled coupling strength between the N M majority 
fermions and the N m = N — Nm minority fermions. We 
will examine in the next sections the following cases 

• a single impurity for the completely imbalanced 
configurations 2:1, 3:1, 4:1 

• a balanced population 2:2 

• a partially imbalanced population 3:2 

the notation being NM'N m . We restrict ourselves to 
these few-body scenarios which are relevant to the ex- 
periment y, Q (the first case is already realized exper- 
imentally [2j) and for which we can obtain reliable ac- 
curate numerical results with the MCTDH method (see 
Appendix). In most cases these are sufficient to reveal 
the mechanisms present also for larger N. 

The special balanced case of two particles (1:1) allows 
for an exact solution not only in ID but also in higher di- 
mensions analyzed in the seminal paper of Busch et. al. 
[34j . The energy curves obtained there solving a tran- 
scendental equation, were also measured experimentally 
in ref. [3] with a very good accuracy via tunneling rates 
after opening the harmonic trap from one side. Extend- 
ing this study to higher particle numbers from the theo- 
retical side is the basic scope of this work. In this context, 
the authors have proposed in ref. [l[ an explicitly cor- 
related pair wave-function (CPWF) for identical bosons 
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which is based on a construction principle inspired by 
the two-body solution (34{ . This CPWF will be extended 
and generalized here to cover the case of two-component 
Fermi mixture. 

The Hamiltonian [Eq. (|T])] containing the harmonic 
trap potential can be separated into a center-of-mass and 
relative motion part. The center of mass R = xi/N 
obeys exclusively in its harmonic oscillator Hamiltonian, 
with the ground state * CM = (iV/7r) 1 / 4 e~ JV ' R2/2 . The 
equation of the relative motion in the general case is much 
more complex, contains kinetic couplings of the relative 
degrees of freedom, includes the interaction terms and 
does not allow for an analytical solution. The CPWF 
which we proposed |l( in the case of bosons for the 
relative motion of N particles consists of products of 
parabolic cylinder functions (PCF) D p (3f|, which is the 
key function occurring in the solution of the two-body 
problem, of the relative coordinate r,j = — Xj of each 
pair. The substantial difference in the case of fermions is 
that the pairs of identical particles (atoms belonging to 
the same hyperfme states) should be put on the levels of 
the Fermi ladder as we will see in the following. 

The generalization of the CPWF we propose here for 
the relative motion reads then: 

P Pm P m 

*c P - n D MnA) n Mprv) n 



*<N M <j 



i<j<N A 



N M <i<j 



(2) 

(up to a normalization constant which we neglect here), 
where P is the number of distinct pairs of atoms be- 
longing to different hyperfme states, Pm and P m are the 
number of pairs of majority and minority atoms respec- 
tively. The pinning of identical fermions belonging to 
the same component on the corresponding Fermi-ladder 
is done by choosing /i = 1 as a parameter for the PCFs 
(see the two last terms of Eq. [5]) which corresponds to 
the energetically lowest orbital of the fermionic relative 
motion: D\(x) = \f2xe~ x / 4 . Therefore the part of the 
total wave function which refers to identical fermions is 
the correct non-interacting ground state resulting from 
the Slatter Determinant: ^ ic i — FJ^ e~ x ^ 2 Yli<j( x i ~ x j) 
with all Xi's and Xj's being exclusively coordinates be- 
longing to the same species of fermions here. For the 
interacting pairs in the wave function the parameter fi of 
the PCF (see first term of Eq. [5]) which depends on the 
interaction strength g should be determined. This can be 
done exploiting the condition 2/3D^ (0) = gDij(0), where 

D'ij = ) and Aj = Ai($ r *i)> which is imposed for 
each pair of atoms belonging to different components in 
the ensemble at r,j = (as we have shown in Ref. [l|) 
due to the discontinuity that is imposed by the delta-like 
interaction potential on first derivative. The resulting 
transcendental equation 



g 2lr(V) 



(3) 



P F(^) 
is solved for /x, selecting the solution in the interval \x € 



[0, 1] which corresponds to the ground state. 

Another key feature of the CPWF is that it reproduces 
the exact analytically known solutions in the limits g = 

and g — > oo, for any N by choosing (3 — \J~^- In these 
limits, /i equals and 1, respectively, and in both cases 



(as well as for every integer fi) D fl (x) 



"He M (a;) 



where He M (x) are the modified Hcrmite polynomials. 
Therefore the total wave function takes the following 
form 
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n »•« n 

i<j<N M N M <i<j 



(4) 



for the non-interacting and 
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Pm 



i<N M <j 



i n ^ n 

i<j<N M N M <i<j 



(5) 

for the infinitely strong repulsion limit [ltl [3(1 HH which 
are the exact solutions in both cases. The CPWF is also 
exact for arbitrary g for the case N = 2 (with f3 = 1) [34| 
and indeed our general Ansatz is inspired by this two 
atom case [H [34[ . 

Using the CPWF as trial function one can treat /3 
as a variational parameter for given values of g and N. 
The form of the pair-correlated function is reminiscent 
of the Bijl-Jastrow form of products, the main differ- 
ence being here that we use the exact solution of the 
two atom case instead of a trigonometric function which 
is the usual approach e.g. in quantum Monte Carlo stud- 
ies. Generalizations of this function can also be extended 
to other cases like bosonic mixtures or even fermion- 
boson mixtures or to many component mixtures. One 
should take \i = 1 for identical fermions and a /i £ [0, 1] 
for interacting pairs according to Eq. [3] Yet one has 
to keep in mind that in case of more than two species 
and for the situation that only a subset of the coupling 
strengths are infinite the above Ansatz may not be ex- 
act. Consequently the cases of a single species bosons 
and two species fermions are the most natural ones for 
which this Ansatz should have its best performance. An- 
other possible generalization is to separate the PCF into 
a product of Hypergeometric and exponential functions 

D fl (x) — 2% e~^U(— ^, g, ^r), and choose two indepen- 
dent variational parameters one in the exponent and one 
in the confluent hypergeometric function of second kind 
U(a,b,x) [35|]. We will restrict ourselves here to the ap- 
proach we addressed above for two-component fermions 
(we will only treat a single case (2:1) variationally with 
a single parameter (3) in order to show the efficiency of 
CPWF as it stands in Eq. [5] to describe the crossover 
from weak to strong interactions. 

A first illustration of the behaviour of the CPWF is 
given for the case 2:1 in Fig. 1 (b) for weak and (c) 
for strong interaction strength, where contour plots of 
the probability density are shown. On the manifolds 
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the identical atoms meet (here the plane x\ = £2) the 
probability density possesses already for weak interac- 
tions [Fig. 1 (b)] a low probability, while stronger in- 
teractions between the non-identical atoms [see Fig. 1 
(c)] reduce the values of the density also on the other 
two planes {x\ = 23, X2 — X3). This is a main intuitive 
picture that the CPWF captures, treating each fermionic 
pair -being interacting or not - correctly at its contact 
point, where the wave function should tend to vanish as 
we increase the coupling strength from weak to strong 
interactions. In the next section the accuracy of this ap- 
proach with respect to the energy is tested by comparison 
with corresponding numerical results for all cases exam- 
ined in this work. The specific behaviour of the densities 
will be addressed in the sections thereafter. 



III. ENERGIES AND ACCURACY OF THE 
ANSATZ 

The energy of the interacting ensemble is a first observ- 
able that can be measured experimentally for example by 
RF-spectroscopy (see e.g. Refs. @, HH, EH) or by deter- 
mining tunneling rates [J] . We compare here the numeri- 
cal MCTDH results to the calculations using the CPWF 
and consequently obtain information on the validity and 
accuracy of the CPWF. 

As discussed in the Appendix MCTDH always involves 
a cut-off for the number of single-particle functions which 
is employed, which in turn imposes also restrictions to 
the number of particles that can be treated with suffi- 
cient accuracy. Most numerical methods are based on a 
truncated Hilbert space, and therefore result in inaccura- 
cies especially when the system is highly correlated and 
a few (even optimal) orbitals are not enough to describe 
it. For the ID case we examine here, there is an efficient 
workaround to overcome this problem, which is related 
to the fact that the energy should converge for g — > 00 
to a certain value equal to the one of the corresponding 
system composed of identical fermions. In Ref. [37[ a 
detailed study of how to overcome the truncated Hilbert 
space problem was performed for a bosonic system in the 
harmonic trap. The main trick is that for a certain set 
of values for the numerical parameters (number of basis 
functions and grid points) one has to rescale g such that 
the energy becomes equal to the fermionic one only at 
g — > 00. Therefore if go is the finite value where the en- 
ergy of the numerical calculation reaches fermionization, 
the corrected rescaled q should be q — -, — ^ — where q n 
is the coupling strength used in the numerical calcula- 
tion. This procedure is available and of importance also 
to our case of two-component fermion system, and we 
use it for extracting the numerically exact energy curves 
shown in Fig. 1 (a) as a function of the inverse inter- 
action strength (analogous to the effective ID scattering 
length, see ref. (18|)- 

In Fig. 1 we also show the energies obtained from 
calculating the expectation value of the Hamiltonian op- 



erator using the CPWF [Eq. fl2J)]. The accuracy of the 
Ansatz is very good especially for weak and very strong 
interactions, yet there are deviations from the numerical 
results in particular for intermediate coupling strengths, 
especially for larger particle numbers (see e.g. the curves 
for the configuration 2:1 in comparison with 4:1). It is to 
be expected that the CPWF works good close to the two 
exactly known limiting cases since it reproduces them 
[see Eqs. ((4]) and The deviation for larger num- 

bers of atoms is due to the two-body nature of the func- 
tion which does not take into account higher correlations. 
Nevertheless a relative error of 2-3% can be achieved also 
by taking variational trial functions with a single varia- 
tional parameter. 

Concerning the physical behaviou, the energy increases 
as the coupling strength gets larger (Fig. 1), but not 
with the same slope for all configurations. The range 
of the energies is bounded by the non-interacting val- 
ues (for the configurations 2:1, 2:2, 3:1, 3:2 and 4:1 it 
is 2.5, 4, 5, 6.5 and 8.5 respectively) and the value at 
fermionization which is Eg^.^ — N 2 /2. This means that 
e.g. the partially imbalanced case 3:2 or the equal pop- 
ulation one 2:2 increase stronger than the corresponding 
impurity cases 4:1 and 3:1 since more interacting pairs of 
atoms are included in the ensemble. On the other hand if 
we consider only the impurity cases, then the gain of en- 
ergy per majority atom (compared to the corresponding 
non-interacting state) at a finite value of the interaction 
strength is larger for a smaller number of atoms. This 
means that by adding majority fermionic atoms the en- 
ergy gain per atom decreases since the majority atoms 
lie higher and higher in energy, such that they almost do 
not feel the impurity fermion at the bottom. Yet when 
the interaction becomes very large the impurity affects 
all atoms and at fermionization the energy gain is for all 
cases equal to Eg^-Eg^ = N 2 / '2- iV|J '2-1/ '2 = N M . 
Further features of the behaviour of the physics of the im- 
purity in the finite Fermi-sea confined in a harmonic trap, 
a case also related to the experiment Q will be illustrated 
in the next section. A relevant question in this context 
is the connection between many-bo dy p olaronic physics 
in the continuum (see eg. Refs. p^ - ll4 l24l [25|) and the 
few-body behaviour of an impurity in the harmonic trap. 



IV. AN IMPURITY ATOM IN A SEA OF 
MAJORITY FERMIONS 

In this section we focus on the impurity physics in ID 
and analyze the corresponding probability densities. Ma- 
jor progress in recent experiments (38j have extended to- 
wards probing the density distribution on a single site 
level. As demonstrated in Figs. 2 and 3 for the config- 
urations 2:1 and 3:1 respectively, the probability density 
shows a vast change for both the majority particles and 
the single minority one. In the same Figures we show 
also the one-body density obtained by the CPWF (Eq. 
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Figure 1: (a) Energy as a function of the inverse interaction strength 1/g for all cases examined in this work (impurity 2:1, 
3:1, 4:1, balanced 2:2, imbalanced 3:2). For each case the lower line corresponds to the numerical results employing MCTDH 
and the line slightly above to the energy computed from the Ansatz CPWF Eq. [2] Contour plots of the density distribution 
X2, 013) | 2 = 0.01 for the configuration 2:1 using the CPWF for a weak (b) g = 0.5 and a strong (c) g — 5.0 interaction 
strength. 




Figure 2: One body densities of (a-d) the impurity and (e-h) majority fermions for the 2:1 case with varying interaction strength 
g using the numerical (MCTDH) approach, and the analytical (CPWF) and variational Ansatz. 
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We focus first on Fig. 2 (a-d) which shows the evolu- 
tion of the one-body density of the impurity from weak 
(g = 0.5,2.5) to stronger interactions (g = 5.0,20.0). 
For weak interactions [g — 0.5 Fig. 2 (a)] the impurity 
density is very close to the ground state orbital of the 
harmonic trap (Gaussian profile). As the interaction in- 
creases [g — 2.5 Fig. 2 (b)] the density becomes broader 
and the peak lowers, which is a typical effect of repul- 
sive interaction also in the case of bosons [l|, [33| . Still 
one may expect that in the case of fermions this single 
impurity atom would be even more confined and pinned 
in the middle of the trap due to the repulsive interac- 
tion with the surrounding majority atoms. However the 
opposite effect is observed: the density is stretched to 
cover a larger space in the trap and becomes delocalised. 
For even higher interaction strengths [g = 5.0 Fig. 2 (c)] 
it even tends to acquire side lobes. More spectacularly 
for very strong interactions, the density of the impurity 
acquires three distinct maxima, an effect which is very ev- 
ident in few-body fermionization profiles also for bosons 
[E [H, Hi| . The number of maxima is equal to the total 
number of atoms in the trap and according to Girardeau's 
theorem [lj| the local properties like densities should be 
exactly the same for a fermionized ensemble and for the 
corresponding ensemble of N identical fermions. This 
vast change of the density of the impurity atom as we 
follow the crossover from weak interactions to fermion- 
ization also for the case of the 3:1 configuration is shown 
in Figs. 3 (a-d). 

The impurity is strongly affected by the repulsion with 
the majority atoms, which in turn are also strongly influ- 
enced by the single minority atom. Their densities for the 
case 2:1 depicted on Fig. 2 (e-h) possess at weak interac- 
tions [g = 0.5 Fig. 2 (e)] two maxima (corresponding to 
a two identical-fermion ground-state profile) which are 
pushed to outwards as the repulsion increases [g = 2.5 
Fig. 2 (f), remember that the corresponding impurity 
density becomes more delocalized in Fig. 2 (a),(b)] and 
start to acquire a plateau and finally an additional peak 
in the middle for very strong coupling [g — 5.0, 20.0 Fig. 
2 (g),(h)]. A similar behaviour is observed also for the 
case 3:1 where the middle peak for weak interactions is 
slightly broadened pushing the outer maxima outwards 
[g = 0.5, 2.5, 5.0 Fig. 3 (e),(f),(g)] and finally is turned to 
two peaks reaching the fermionization profile for strong 
coupling [g = 20.0 Fig. 3 (h)]. The impurity fermion af- 
fects therefore the density profile of the majority atoms 
equally strong, leading to the extreme infinitely repulsive 
limit where it becomes effectively indistinguishable from 
the majority fermions (strictly only for the local proper- 
ties, of course). 

An issue that is worth noting here is that the very 
strong change in the density profile (where more distinct 
peaks appear) is happening rather abruptly for a very 
high interaction strength. The latter distinguishes the 
case of few fermions from that of few bosons which ac- 
quire substantial peaks already for intermediate interac- 
tions (see e.g. [l|). In other words, this is a characteristic 



behaviour of two-composite few fermion systems, which 
becomes interesting in the very strong correlation regime, 
being the important regime of ID physics. This regime 
is difficult to access both for the experiment (one has to 
approach very closely to the g — > oo resonance, and si- 
multaneously avoid losses to the molecular side |2j) and 
for numerical calculations since a lot of basis functions 
have to be included, rendering the treatment of a result- 
ing huge Hilbert space very difficult. The application 
and convergence of MCTDH for such strong interactions 
is a very subtle issue, and has been tested extensively 
in this work, leading to the reliable results for the cases 
presented, which we also compare with the correspond- 
ing CPWF (which naturally provides the correct profile 
for infinite repulsion, thus not running into the affore- 
mentioned difficulties). We note that here the values of 
g are the actual ones used in the calculation and are not 
rescaled like in the case of the energy. 

The accuracy of CPWF to describe the evolution of the 
densities of the impurity or the majority particles is very 
good close to the two, the non-interacting and fermion- 
ization, limits as expected and describes the crossover 
qualitatively very well. Yet at intermediate interactions, 
the CPWF densities tend to be broader and acquire mi- 
nor additional peaks compared to the numerical results 
[see Figs. 2, 3 (b),(c),(f),(g)]. This is in accordance 
with behaviour of the corresponding energies, i.e., the 
fact that the analytical CPWF energies lie above the nu- 
merical ones in the intermediate interaction regime [see 
Fig. 1(a)]. For the 2:1 case we present also the densities 
obtained by a variational treatment of the parameter (3 
(minimization of the energy using the CPWF ^ cp as a 
trial function). We observe in Fig. 2 that the variational 
treatment of the CPWF represents a significant improve- 
ment over the fixed j3 CPWF case: the corresponding 
densities get closer to the numerical ones. Further im- 
provements of the Ansatz to capture the exact quantita- 
tive behaviour may include more variational parameters 
or a different construction in terms of hypergeometric 
functions discussed in section II. Yet we may emphasize 
that the correct qualitative behaviour is captured by the 
explicit function ^ cp . A similar good agreement exists 
for the next two cases of partial and complete balance of 
populations that we present in the next section. 



V. EQUAL POPULATIONS AND PARTIAL 
IMBALANCE 

In the following we analyze the density profiles of two 
cases which exhibit a different crossover to fermioniza- 
tion than the impurity physics discussed above. Since 
the Ansatz of the CPWF possesses a similar accuracy as 
already depicted in Fig. 1 (a) we will concentrate here 
on the physical effects for these cases and present only 
the numerically obtained densities in Fig. 4. 

Fig. 4 (a) shows the one-body density profiles for dif- 
ferent interaction strengths for the case of equal popula- 
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tions (2:2). Note that the two hyperhne states possess 
exactly the same profile features in this completely bal- 
anced case. The two initial peaks for weak g = 0.5 get 
broadened and flattened (g — 2.0) and start already for 
intermediate g = 5.0 to acquire side peaks, resulting in 4 
distinct peaks for a strong coupling g — 15 according to 
the fermionization theorem. The crossover is here much 
smoother than in the case of the impurity (reminiscent of 
the bosonic case [H, HH) with each peak (corresponding 
to a pair) spiting to two. The behaviour for equal popula- 
tions can be also generalized to higher particle numbers, 
and we expect no substantial differences which we have 
checked for the few body case of 3:3 where the 3 initial 
peaks evolve to 6 again due to repulsive pair splitting. 

Another very interesting scenario is that of a partially 
imbalanced population, which we illustrate here for the 
case 3:2 in Fig. 4 (b) and (c) showing the density pro- 
files with increasing coupling strengths for the minority 
and the majority atoms, respectively. The two minor- 
ity atoms possess 2 lobes for weak interactions g = 0.5 
which are slightly pushed outwards in the trap for inter- 
mediate interactions g — 5.0, 10 forming a small plateau 
in the middle. Only for a strong coupling g — 14 a third 
maxima rises in the middle of the potential and some 
outermost side wings are forming, which result to two 
additional peaks at g = 16. We see again here a rather 
abrupt and interesting behaviour for very strong inter- 
actions, like in the case of the impurity fermion. The 
majority atoms [Fig. 4 (c)] also hold the three-lobe pro- 
file for intermediate interactions (g = 0.5,5.0,10), while 
from the plateaus forming between the peaks two addi- 
tional peaks emerge near g — 14.0 resulting in a profile 
with five maxima. Inspecting Figs. 4 (b) and (c) the two 
peaks of the minority density (g = 10) which are located 
at the positions of the two plateaus of the majority den- 
sity seem to represent a rather stable configuration for 
the repulsive particles. Still the tendency to fcrmionizc 
in order to avoid each other as much as possible breaks 
this configuration giving rise to the additional peaks. 

The above mechanisms of the crossover to fermioniza- 
tion are present also in ensembles with larger number of 
atoms, yet they are more prominent and illustrative for 
few-atom systems. The connection between these mech- 
anisms and the many-body phases like BCS and FFLO, 
is an open research pathway. 



VI. CONCLUDING REMARKS AND OUTLOOK 

We have employed exact numerical calculations 
(MCTDH) and an analytical approach (correlated-pair 
wave function CPWF) to explore the properties of few- 
fermion two-component mixtures in an one-dimensional 
harmonic trap. Our Ansatz for the function of the rela- 
tive motion CPWF which is here generalized to cover the 
case of mixtures, is shown to be in a very good agreement 
(for the energies and the densities) with the numerical 
results of MCTDH (used here for the first time for the 
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Figure 4: one body densities for the case (a) 2:2 (the densities 
for the two components are identical) and for the case 3:2 for 
the (b) minority and (c) majority fermions. 
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few-fermion case). We examined three general cases: (a) 
an impurity in a sea of majority fermions, where we have 
shown a strong variation of the density with an increasing 
strength of the repulsive interaction both for the impurity 
atom and the majority ones, (b) an equal population case 
for which the behaviour is smoother with two peaks aris- 
ing in the place of one as the coupling strength increases 
for each pair of atoms and (c) a partially imbalanced case 
where very close to the resonance additional peaks arise 
on top of the plateaus in the density. 

Most of these theoretical considerations are also stud- 
ied in a very controllable manner in experiments of few 
fcrmion systems with a deterministically preparable atom 
number [2|-l4| . The generalization of our Ansatz for the 
CPWF to mixtures is open to further improvements and 
applications which we initiated and sketched in this work, 
as well as to the extraction of other observables unrav- 
eling mechanisms and physical properties characterizing 
few-body systems in connection with and analogue to the 
corresponding many-body physics. 
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Appendix A: Computational Method 

In this work we have used for the numerical calcula- 
tions the Multi- Configurational Time-Dependent Hartree 
(MCTDH) method [Hi, H§ E3| , which is a wave-packet 
dynamical tool with an outstanding efficiency for high- 
dimensional systems. The method has been applied suc- 
cessfully to study few-boson dynamics and their sta- 
tionary properties (33{. Here we apply this method for 
first time to fermionic systems. Later developments of 
MCTDH include by construction the bosonic (MCT- 
DHB) 53 or the fermionic (MCTDHF) {44] character 
of the system. We resort here to the initial Heidelberg 
MCTDH package [41| implemented for distinguishable 
particles and enforce the correct antisymmetry by the 
corresponding choice of the coefficients of the MCTDH 
expansion. 

The underlying idea of MCTDH is to treat the time- 
dependent Schrodinger equation 



*(Q,0) = * (<2) 



(Al) 



Hartree products $j: 



Mi Tlf 



K=l 



where J = (ji . . .jf) are the different configurations and 
f = N denotes the number of degrees of freedom and 
Q = (x\, . . . , Xf) T . The single-particle functions SPFs 

ipj^ are in turn represented in a fixed, primitive basis 
implemented on a grid. For indistinguishable particles 
as in our case within the same component of fermions, 
the single-particle functions for each degree of freedom 
are chosen to be the first Nm or N m harmonic oscilla- 
tor states (identical for each atom belonging to the same 
component). 

In the above expansion, both the coefficients Aj and 
the Hartree products <I>j are time-dependent. Using the 
Dirac-Frenkel variational principle, one can derive equa- 
tions of motion for both Aj, <I>j. This offers a huge ad- 
vantage since the basis {$,/(■, t)} is variationally optimal 
at each time t, and thus we can choose a small number of 
SPFs. The correct permutation symmetry which is im- 
portant here can be enforced on the coefficients Aj such 
that the total expansion is antisymmetric with respect to 
permutation of two atoms of the same species. 

In addition the Heidelberg MCTDH package [4l[ incor- 
porates the so-called relaxation method which provides a 
way to obtain the lowest eigenstates of the system by 
propagating some wave function by the non-unitary 
e~ Hr propagation in imaginary time. As t -> oo, this 
automatically damps out any contribution but that stem- 
ming from the true ground state, 



-Ejt 



IWI*0>. 



as an initial-value problem by an expansion in terms of 



In our case it is important that the ground state is not 
as in the bosonic case given for free as the lowest state of 
indistinguishable particles system, but we have to apply 
a more sophisticated scheme termed improved relaxation 
[42J where one starts for g — with the two-component 
fermionic ground state placing the fermions on the Fermi 
ladder correctly. Here (^f\H — E\^) is minimised with 
respect to both the coefficients Aj and the configurations 
$j. The equations of motion are solved iteratively, first 
for A j {t) (by diagonalisation of ($j\H\$k) with fixed 
3>j) and then propagating <i>j in imaginary time over 
a short period. The cycle will then be repeated until 
convergence is achieved. 

The computational effort of this method scales expo- 
nentially with the number of degrees of freedom, n N . 
Just as an illustration, using 15 orbitals and N = 5 
requires 7.6 • 10 5 configurations J. This restricts our 
analysis in the current setup to 5 atoms totally, given 
that we need to cover the whole range of interactions 
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and especially the stronger coupling regime. By con- 
trast, the dependence on the primitive basis, and thus 
on the grid points, is not as severe. In our case, the grid 
spacing should of course be small enough to sample the 
interaction potential. Since the truncation of the Hilbert 
space (to 15 or 20 orbitals for N — 4,5 respectively) is 



unavoidable one may employ additional tricks that heal 
this problem, as the rescaling of the interaction param- 
eter proposed in [13] . The ID systems with the known 
fcrmionization limit allow for such a treatment that we 
employ here as explained in more detail in Section III. 
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